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ABSTRACT 

A method is described for calculation of the radiation pattern of large aperture 
antennas. A piece-wise linear approximation of the aperture field using overlapping 
pyramidal basis functions allows the radiation pattern of an aperture antenna to be 
calculated as though it were a two-dimensional array. The calculation of radiation 
pattern data versus 0 and <p, suitable for 3-D or contour plot algorithms, is achieved 
by locating the array in the yz-plane and performing a summation over the aperture 
field data sampled on a square grid. A FORTRAN subroutine is provided for 
performing radiation pattern calculations. Numerical results are included to 
demonstrate the accuracy and convergence of the method. These numerical results 
indicate that typical accuracies of ±0.idB for Directivity, ±ldB for the 1st Sidelobe 
Level, and ±2dB for the 2nd Sidelobe Level can be obtained with an aperture grid of 
45x^5 points and requires approximately 0.02 seconds CPU time per far-field data 
point on a VAX 1 1/750 with a floating point accelerator. 
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INTRODUCTION 


This paper describes a method for calculating the radiation pattern of a large 
aperture antenna. In order to avoid the necessity of performing a double numerical 
integration to calculate the radiation patterns, overlapping pyramidal functions are 
used to develop a piece-wise linear representation of the aperture field. This 
representation allows the radiation pattern of the aperture antenna to be calculated 
as though the antenna is a two-dimensional array antenna with an array element 
pattern described by the radiation pattern of a small square aperture with a 
pyramidal distribution. The problem of calculating the radiation pattern of the large 
aperture then becomes one of calculating the array factor of a two-dimensional array 
and multiplying by the array element pattern. The advantage of this approach is that 
the double numerical integration becomes a summation over the complex amplitudes 
of the aperture field at the sampled grid points weighted with an exponential factor 
corresponding to the grid point location. The problem is formulated such that the 
array (i.e., aperture) is located in the yz-plane. When the aperture is in the yz- 
plane, the far-field radiation pattern calculation over 0<J>-space can be accomplished 
more efficiently. 

The present problem is formulated by first expressing the far-field radiation 
pattern of an aperture located in the yz-plane in terms of the Fourier transforms of 
the tangential aperture electric fields. A piece-wise linear representation of the 
aperture field distribution is obtained by the superposition of 3 two-dimensional 
arrays of pyramids. The three arrays are interleaved such that the superposition 
yields a 3-dimensional surface approximated by triangular facets. The Fourier 
transform of this distribution is expressed in a series form which allows an 
efficient computer code to be written for the radiation pattern calculation. In order 
to implement this approach, the aperture field must be sampled over a rectangular 
area with the same sample spacing in both directions. If the aperture boundary is not 
rectangular, the field values at the sample points outside the aperture must be set to 
zero. 
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SYMBOLS 


VyV 

B m«yy 

E a (y,d 


E & E 

y z 
E e & 


nm 

E A &e b 

nm nm 

f(k ,k ) 
y’ z 

f A (k ,k ) 

o y’ z 

f B (k ,k ) 

o y’ z 

f (k ,k ) 
y y’ z 

f (k ,k ) 
z l y’ z 


F e &F * 


Grid spacing for sampled aperture field data. 

Series defined in equations (28a) and (29a). 

Series defined in equations (28b) and (29b). 

Electric field distribution in aperture. 

Vector electric field in aperture, 
y- and z-components of E q . 

0- and ^-components of electric Field in far-field region. 

Complex amplitude of aperture electric field at sample point (y^,z°) 
Definition of type-A and type-B pyramidal basis functions. 

Fourier transform of aperture electric field distribution. 

Fourier transforms of type-A and type-B basis functions. 

Fourier transforms of y- and z-components of aperture electric field. 

Electric vector potential function. 

0- and (^-components of F. 

=\fT 


k 

k &k 

y z 

m 

M 


=2n/X, Wave number in free space. 

Fourier transform variables in y- and z-directions. 
Grid index in y-direction. 

Total number of sample grid lines in y-direction. 
Equivalent vector magnetic current in aperture. 
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Grid index in z-direction. 
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n 

n 

N 

r,6,<t> 

r o 

S 

a 
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x.y.z 

y & z 



o 

z 

n 


X 


Unit vector normal to aperture. 

Total number of sample grid lines in z-direction. 
Coordinate variables in spherical coordinate system. 
Radius of circular aperture. 

Area of aperture. 

Width of square aperture. 

Coordinate variables in cartesian coordinate system. 

Unit vectors in y- and z-directions. 

Variables defined in equation (14). 

Variables defined in equation (15). 

y-coordinate of m-th sampling grid line in aperture plane, 
z-coordinate of n-th sampling grid line in aperture plane. 
Wavelength in free space. 
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METHOD 


The electric field intensity in the far-field region of an antenna is given by 

E 0 = 'J k F <p 

E d> = J k F 0 & 

where Fq and are the spherical components of the vector potential function 


F = 



(3) 


where the equivalent magnetic current, M g , is obtained from the tangential electric 
field, E , in the aperture (x=0) 

M s = -2n x E a (4) 

M s = 2 (E z y - E y z) (5) 


By substituting 

ky = k sin0 sin$ 

k = k cos0 
z 


the spherical components of the radiated field become 


E 0 = ZJ ^F~ t ? z {k y' k z ] cos<t> l 

E (f> = ^ r " t f Z ( k y’ k Z 5 COS0 Sin< ^ + f y( k y> k z ) Sin0 ] 


( 6 ) 

(7) 

( 8 ) 
(9) 


where f^(k ,k z ) and f z ( k y 5 k z ) are the Fourier transforms of the aperture electric 
fields defined as 


rr jk y jk z 

y k y> k z) = JJ E y (y,z) e y e z dy dz 


( 10 ) 
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(il) 


f z (k y ,y = If E z (y, Z ) .-"y 


j V jk z z , a 

e 1 e dy dz 


In order to calculate the radiation fields, the integrals in (10) and (11) must be 
evaluated. For practical aperture field distributions, these integrals have to be 
evaluated numerically. In the case of sampled aperture data, in which the data are 
on a rectangular grid, Fast Fourier Transforms can be used to evaluate the far-field. 
However, use of FFT algorithms results in far-field data which are also spaced on a 
rectangular grid in the transform variables, k^ and k^, instead of 0 and <p as is 
sometimes desired in order to utilize existing 3-D or contour plotting algorithms. 

Another approach would be to expand the aperture field in a set of basis functions 
whose Fourier transforms can be evaluated in closed form. The basis function used 
here is a pyramid with a square base, a height equal to the sampled aperture field 
and with the apex located at the sampled aperture data point. Use of the pyramidal 
pulse basis function allows the aperture field to be approximated by a piece-wise 
linear representation in three-dimensions. In order to accomplish this piece-wise 
linear representation, two types of pyramidal basis functions are needed as defined by 

type-A : 


E A = 

nm 


o 

nm 

o 

nm 

o 

nm 

o 

nm 


type-B : 


E B = 
nm 


r° 

'nr 

: 0 
'nm 

7 ° 

'nm 

:° 

'nm 


(1-V=) 

(O^z ^a ; -z ^y ^z ) 


{O^y ^a ; -y ^z ^y } 
1 y m n 7 m J 

( 1 + z„/a 1 

{-a^z ^0 ; z £y £-z } 
1 n ’ n y m n J 

(1 *y^ /e) 

l-aSy m *o i y m %s-y m ) 

1 1 • ( ym +z n )/a ] 

(°Sy m <: a ; 0Sz n S(a-y m )) 

[ 1 - (y - z )/a ] 
1 k/ m n' 1 

(° S y m Sa ; - (a '>'m liz n S0 l 

i 1 + <y m +z ^ /a ) 

!-^y m *0 ; -(a+y m )Sz n S0) 

i 1 + <V z n )/a i 

(- a Sy m i° ; 0Sz n S(a+y ni )} 


( 12 ) 


(13) 
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where, 


y m = y-y 


m 


z = z - z„ 

n n 


d4) 

(15) 


The location of the sampled aperture data point is defined by (y^> z ^)» the amplitude 
of the sampled aperture field is represented by the complex variable E^, and a 
represents the sampling interval in the aperture plane. 

Figure 1 demonstrates the use of these pyramidal basis functions to approximate a 
cos(7ry/w)cos(7rz/w) distribution using 25 (i.e. 5x5) sample points. The piece-wise 
linear approximation is obtained by the superposition of two interleaved arrays of 
type-A basis functions with one array of type-B basis functions. This interleaving of 
arrays of pyramidal basis functions results in an aperture distribution given by 


e ( y ,z) = 2 2 e a + 2 2 e a + 22 e b + 2 2 e b 


m n 
0 0 


m n 
e e 


nm 


m n 
o e 


nm 


m n 
e o 


nm 


(16) 


where m and n indicate even values of m and n; whereas, rn and n indicate odd 
e e ’’oo 

values of m and n. The Fourier transform of (16) is given by 


jk y jk z 

f(>V k z ) = 2 2 E nm ■ // d'V 3 ) e y e d ydz 

m " jk y J k z 

+J/(i-y rri /a) e y e z dy dz 


m n 
o o 


+// (1+z /a) e y ' 


j k v y K z 


dy dz 


jk y jk z -v 
+// (l+y m /a) e y e z dy dz | 


r jk v jk z 

+ 2 2 E ^n { // d-V 3 ) e y e z d y dz 

m n V . 


m n 
e e 


j k y jk z 

+// (l-y m /a) e y e z dy dz 


j k v y JM 


+//( i+z /a) e y e z dy dz 


+JJ d+y m /a) e y ' 


J k v y JM 


dy dz | 
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+ 2 2 ^ { // [ 1 -(y m +z n )/al * yY z dy dz 

m o n e jk y jk z 

+//( 1 -(y m ' z n )/ a ] e y e z d y dz 

jk y jk z 

+//[ 1 +(ym +z n )/a l e y e d y dz 

jk y jk z -v 

+//[i+(y m -z n )/^ e y e z d y dz ] 

+ 2 2 { // [l-(y m + z n )/a] e j ^ Z dy dz 

m e n o jk y jk z 

+//[Hy m -z n )/a] e y e z d y dz 

jk y jk z 

+//[ 1 +(y r n +z n )/^ e y e z d y dz 

jk y jk z -v 
+//[i + ( y m -z n )/ a ] e y e z dy dz j 

Using (14) and (15), equation (17) becomes 

o o 
d v y m d z z n 



A R 

where f (k ,k ) and f (k ,k ) are defined as 
o y’ z o y’ z 


f A (k ,k ) = 
o y’ z' 


a z 
r r n 


M j^y^m ^z z n , , 

( 1 V a)e e 


0 -z 


ik y ik z 

(1-y /a) e y m e z n dz dy 
7 m n 7 m 


0 -y 


m 
0 -z 


+ 


n ik y ik z 

(i+z n /a) e y m e Z n dy m dz n 


-a z 

n 

0 -y 

r r } ] 


m 


MJ . V" ^ k z Z n j . 
(i+y m /a) e 1 e dz dy, 


n y m 


- a y 


m 


(19) 


f B (k ,k ) = 
o y’ z 


a (a-y ) 

r r 


r ' m ik y ik z 

! t 1( y m +z n )/a 1 e yme " dz n d ^m 


0 0 
a 0 

* 

•( a -y m ) 
0 0 


ik y jk z 

[ l-(y -z )/a ] e y m e z n dz dy 
1 7 m n 1 n 7 m 


-a -(a+y m ) 

P , (a+y m ) 


[ l+(y +z )/a ] e y m e z n dz dy 
1 v m n 1 n 7 m 


-a 


[ l+(y m -z n )/a ] e yym e^ z n dz R dy : 


0 


( 20 ) 


f A (k ,k ) and f B (k ,k ) are the Fourier transforms of (12) and (13) with E° =1 
o y z o y’ z nm 
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and y°=z°=0. After much algebraic manipulation, f A (k ,k ) and f B (k ,k ) can be 
y m n ° r ’ o y’ z o y’ z 

expressed as 


f o = 


4a" 


i (V) 2 - M 2 J 


sin(k a) cos(k z a) sin(k z a) cos(k a) 

Vi Fi Z ~ 


( 21 ) 


f B (k ,k ) 

o ' y’ z 


4a" 


l (V) 2 - M 2 


sin(k z a) 

k a 
z 


sin(k a) 

fi" 


( 22 ) 


with limiting cases given by 

2 

f 0 VV = 2 (H ( 1 - sin(2k y a) /(2 k y a) ) 


rB 


(y 

4 

\ 

fi 1 

2 f 

k 


t y J 

V 


(23) 


(24) 


(0,0) = 


4a" 


f® (0,0) = 


(25) 

(26) 


The Fourier transform of the aperture field can now be written as 
M 


f(k 


• ° i 

1V1 r a n 'l jy k 

= 2 =1 [ A m (k y ,k z ) f o A (k y ,k z ) + B m (k y ,k z ) f o B (k y ,k z ) j e m y (27) 
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where, for m=l,3,5,...,M 


Wz> = 

N jz°k 

2 E° e nz 
^ nm 
n=l 

(n odd) 

(28a) 

B m »y k z > = 

N-l jz°k 

c o J n z 
> E e 
nm 

n=2 

(n even) 

(28b) 

whereas, for m=2,4,6 

, ... ,M-1 



A mVz> = 

N-l jz°k 

r-0 J n z 
> E e 
nm 

n=2 

(n even) 

(29a) 

B m Vz> = 

N jz°k 

■ST r-o J n z 
> E e 
■^4 nm 
n=l 

(n odd) 

(29b) 


and where both M and N are assumed to be odd. 

Using (27) to represent the aperture field transforms in (8) and (9) allows the 
far-field pattern of the large aperture to be calculated from interleaved arrays, 
whose element patterns are obtained from (21) and (22). 



RESULTS 


A FORTRAN subroutine for calculating the far-field pattern of an aperture is 
given in the Appendix. This subroutine was used to evaluate the accuracy and 
convergence of the method by calculating the far-field patterns for circular aperture 
field distributions specified at discrete points on a square grid of NxN grid lines. 

The distributions used in the evaluation were obtained by sampling the expression 
Ey(r)=C+(l-C)[i-(r/r 0 ) 2 ] over a square grid for values of C equal to 0.0, 0.316 and 
1.0, as shown in figure 2. These piece-wise linear aperture distributions are shown 
in figures 3-5 for grid densities of 15x15, 35x35 and 55x55. Three-dimensional 
plots of the radiation patterns for an aperture diameter of 1000 A are shown in 
figures 6-8. These patterns exhibit the characteristic ring sidelobes expected from 
a circular aperture with a r-dependent excitation distribution; however, careful 
examination of the patterns for N=15 shows that each of the sidelobes, after the 
first, has an undulating character as it encircles the main beam. This becomes 
more obvious for the contour radiation patterns in figures 9-11. Since very little 
change is noticed between N=35 and N=55, these data tend to indicate that a grid 
density of 35x35 may be sufficiently dense for accurate radiation pattern 
calculations; however, one should be cautioned about the possible occurrence of 
"grating" lobes in the calculated patterns caused by large spacings between aperture 
sample points. The angular regions where these grating lobes can occur are 
predictable and should be avoided when calculating radiation patterns. 

In order to better establish the minimum grid density needed, the radiation 
pattern parameters of directivity, beamwidth and the first two sidelobe levels were 
calculated in the 0=90° plane over a wide range of grid densities. These data are 
plotted in figures 12-15. These data indicate that very accurate radiation patterns 
can be calculated with this method by using an aperture grid density of only 45x45 
sample points; however, one must be reminded that this convergence criterion was 
established for the class of circular aperture distributions defined as a parabola-on-a- 
pedestal and is intended to only be a guide as to the grid density necessary. One 
should always check the convergence of the sampling density for the type of 
distribution one expects for his particular application. 


12 



In order to determine the computational efficiency of the method, far-field 
patterns were calculated over a square grid in 0 (^-coordinates (81x81=6561 far- 
field points). The average CPU time per far-field point is plotted in figure 16 
versus the number of aperture grid lines. The CPU time in figure 16 was obtained 
on a VAX 11/750 with a floating-point accelerator installed. The corresponding 
CPU time on an 1BM-PC with an 8087 installed is plotted in figure 17. It was 
observed that the average CPU time per far-field point remained nearly constant as 
the number of far-field points increased. A note of caution: the CPU time per far- 
field point will probably be different for a 00 grid which is not square (i.e. the 
number of 0 and 0 pattern angles are not equal). The CPU times in figures 16 and 
17 were each based on a single subroutine call and includes the time required to 
calculate both the 0- and 0-components of the far-field from both the y- and z- 
components of the aperture field. 
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CONCLUSION 


An efficient numerical method of calculating the volumetric (i.e. over a 00- 
region) radiation pattern of large aperture antennas has been described. Examples 
have been given which demonstrate that very accurate results can be expected with 
an aperture field sample grid of 45x45 points and requires only 0.02 seconds CPU 
time per far-field point on a VAX 1 1/750 with a floating-point accelerator installed. 
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Interleaved arrays of pyramidal basis functions to 
form piece-wise linear representation of 3-D surface 
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Figure 2. Circular aperture -field distributions. 
E y (r) = C + < 1-C) [ l-<r/r 0 )2] 
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Figure 4. Piece-wise linear representation of aperture field. 
Ey<r> = 0.316 + U-0 .316) [ l-<r/r 0 )2 
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Figure 5. Piece-wise linear represen tat i on of aperture field 
Ey(r) = 1 
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Figure A. Three-dimensional plots of calculated far— field 
patterns for aperture distributions of figure 3. 
< 2r 0 =l O00X) 
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■> Figure 7. Three-di mens i onal plots of calculated fai — field 

patterns for aperture distributions of figure 4. 
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Figure 9. Contour plots of calculated far— field patterns 
for aperture distributions of figure 3. 

(5dB increments) (2r o -1000X) 
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Figure 10. Contour plots of calculated far-field patterns 
for aperture distributions of figure 4. 

<5dB increments) <2r o =100QX> 
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Figure 11. Contour plots o-f calculated far-field patterns 
-for aperture distributions o-f -figure 5. 

<5dB increments) <2r«=1000X.) 







Figure 12. Calculated far-field directivity versus 
circular aperture grid density <M=N> . 
Ey<r )=C+ < 1-C) [ 1 -<r/r 0 )2] 
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Figure 13. Calculated far-field 3dB beamwidth versus 
circular aperture grid density <M=N) . 
E y <r)=C+<l-C)[ l-(r/r 0 )2] 
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Figure 14. Calculated far-field 1st sidelobe level versus 
circular aperture grid density <M=N> . 

Ey<r >=C+< 1 -C) C l-<r/r 0 )2] 







Figure 16. VAX 11/750 CPU time for 81X81 far-field point* 
versus aperture grid density (M=N> . 
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APPENDIX 


« 


Subroutine RADPAT 

This subroutine calculates the radiation pattern over a specified Q<p angular region 
in the far-field of an aperture. The amplitude and phase of the aperture electric 
field distribution are specified at discrete points over a rectangular area which 
encloses the aperture boundaries. The discrete field values in the aperture plane are 
constrained to have the same sample spacing in both directions and the number of 
samples in both directions must be an odd number. The parameters in the CALL 
statement to the subroutine are defined by comment statements at the beginning of 
the subroutine listing. 


* 

( 
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oooooooooooooonoooooooooooooooooooo 


SUBROUTINE RADPAT < DD ,NN,MM,NZ ,MY, EYM , EYP , EZM , EZP , AA , 

+ T1 ,T2,T3,P1 f P2,P3) 

C*************************************************************** 
COMPUTES RADIATION PATTERN FOR APERTURE WITH DISTRIBUTION * 
SPECIFIED AS POINTS ON SQUARE GRID IN YZ-PLANE. * 

* 

METHOD AND PROGRAM BY: M. C. BAILEY * 

* 

THETA, PHI, DBT AND DBP OUTPUT TO UNIT #12 * 

* 

STEP ANGLE = THETA * 

SCAN ANGLE = PHI * 

BORESITE DIRECTION IS (THETA=90 ,PHI=0) * 

ft*********************************************************** 

DD * GRID SPACING IN WAVELENGTHS. < SQUARE GRID) * 

NN = NUMBER <ODD> OF APERTURE FIELD POINTS IN Z-DI RECTI ON* 

MM * NUMBER (ODD) OF APERTURE FIELD POINTS IN Y-DI RECTI ON* 

NZ ,MY = ARRAY DIMENSIONS (SEE NOTE BELOW). * 

* 

APERTURE ELECTRIC FIELD AT ( Z-SUBN,Y-SUBM) : * 

EYM(N,M) = Y-COMPONENT MAGNITUDE * 

EYP(N,M> = Y-COMPONENT PHASE (RADIANS) * 

EZM(N,M> = Z -COMPONENT MAGNITUDE * 

EZP(N,M> = Z -COMPONENT PHASE (RADIANS) * 

* 

AA * TEMPORARY STORAGE ARRAY. * 

* 

PATTERN ANGLES (DEGREES): * 

T1,T2,T3 = THETA ( INITIAL , FINAL, INCREMENT) * 

PI ,P2,P3 = PHI (INITIAL, FINAL, INCREMENT) * 

* 

NOTE: IN CALLING PROGRAM; * 

EYM, EYP, EZM, EZP, AA * 

MUST BE DIMENSIONED AS FOLLOWS: * 

* 

EYM(NZ.MY) ,EYP(NZ,MY) * 

EZM(NZ ,MY) ,EZP(NZ ,MY> * 

AA(8,MY> * 

C**************************************************************» 
DIMENSION EYM(NZ ,MY> ,EYP(NZ ,MY) 

DIMENSION EZM(NZ ,MY) ,EZP(NZ ,MY) 

DIMENSION AA( 8 ,MY> 

PI=2.*ASIN(1 .0) 

PI 02=0 . 5*PI 
DT0R=PI/1 80 . 

AK«2.*PI*DD 

MM1=MM-1 

NN1=NN-1 
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IF< <MM1 -2*<MMl/2> > . EQ.O> GO TO 20 

MM-MM+1 

MM1-MM-1 

IF<MM.GT.MY> GO TO 300 * 

DO 10 N=1 ,NN 
EYM(N,MM)=0.0 
EYP<N,MM>=0.0 
EZM<N,MM)=0 .0 
EZP<N,MM>*0.0 
10 CONTINUE 

20 IF< <NN1 -2*<NNl/2> ) . EQ.O) GO TO 40 

NN=NN+ 1 
NN1==NN-1 

IF(NN.GT.NZ) GO TO 300 
DO 30 M=1 ,MM 
EYM<NN,M)»0.0 
EYP<NN,M>“0 .0 
EZM<NN,M>»0.0 
EZP(NN,M)=0 .0 
30 CONTINUE 

40 NT«<NN+l)/2 

MT=<MM+ 1 )/2 
DDSQ=DD*DD 
C0=4.0*DDSG> 

C2=2.0*DDSQ/3.0 

C1=2.0*C2 

UO-1 .0 

PA=0.0 

PB=0.0 

DO 70 N=1 ,NN,2 
DO 50 M*1 ,MM ,2 

PA=PA+ < EYM<N , M) *EYM<N ,M) +EZM(N ,M> *EZM<N ,M> ) 

50 CONTINUE 

DO 60 M=2,MM1,2 

PB»PB+<EYM<N,M>*EYM<N,M>+EZM<N,M>*EZM<N,M>> 

60 CONTINUE 

70 CONTINUE 

DO 100 N=2,NN1,2 
DO 80 M=*2,MM1,2 

PA=PA+ < EYM(N ,M> *EYM<N ,M> +EZM(N ,M> *EZM<N ,M> ) 

80 CONTINUE 

DO 90 M=1 ,MM,2 

PB=PB+<EYM<N,M>#EYM<N,M>+EZM<N,M>*EZM<N,M>> r 

90 CONTINUE 

100 CONTINUE 

PRD“< Cl *PA+C2#PB> < 

110 U0=PRD/<4.*PI> 
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DO 120 M=1 ,MM 
DO 120 N=1 jNN 
A=EYM<N,M> 

P=EYP(N,M) 

EYM<N,M>=A*COS<P> 

EYP<N,M>=A*SIN<P> 

A=EZM<N,M> 

P=EZP<N,M> 

EZM<N,M>=A*COS<P> 

EZP<N,M>«A*SIN<P> 

120 CONTINUE 

NUMTH=<T2~T1 )/T3+ 1 . 5 
NUMPH=(P2-P1)/P3+1 .5 
WRITE < 12,130 >NUMTH ,NUMPH 
130 F0RMAT<2I5> 

DO 290 I»1 ,NUMTH 
THETA=Tl+< 1-1 >*T3 
TH=THETA*DTOR 
CT=COS(TH) 

ST=SIN<TH> 

AKZ=AK*CT 
AKS=AK*ST 
SAZ»1 .0 

I F<ABS<AI<Z ) . GT. 1 .E-05)SAZ ss SIN<AKZ)/AKZ 
CAKZ=COS< AKZ ) 

AKZSQ*=AK2*AKZ 
DO 140 M=1 ,MM 
DO 140 N*1 ,8 
AA(N,M)=0.0 
140 CONTINUE 

DO 170 N=1 ,NN,2 

EN=<N-NT>*AKZ 

ER=»COS(EN) 

EI=SIN< EN) 

DO 150 M=1 ,MM,2 

AA< 1 ,M>=AA<1 ,M>+<ER*EYM<N,M>-EI*EYP<N,M>> 
AA< 2 ,M>=AA< 2 ,M> + < El *EYM<N ,M) +ER*EYP<N ,M> ) 
AA< 3 f M>=AA<3,M) + <ER*EZM<N,M>-EI *EZP<N,M> ) 
AA( 4 ,M>=AA< 4 ,M> + < El *EZM<N ,M> +ER*EZP<N ,M> ) 
150 CONTINUE 

DO 160 M=2,MM1 ,2 

AA< 5 ,M>=AA< 5 ,M> + < ER*EYM<N ,M> -El *EYP<N ,M) ) 
AA<6,M>=AA<6,M)+<EI*EYM<N,M>+ER*EYP<N,M>> 
AA<7,M>=AA<7,M>+<ER*EZM<N,M>-EI#EZP<N,M>> 
AA<8 f M>=AA<8,M>+<EI*EZM<N,M>+ER*EZP<N,M>> 
160 CONTINUE 
170 CONTINUE 
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DO 200 N«2,NN1 ,2 

EN=<N-NT)#AKZ 

ER=COS<EN> 

EI=SIN(EN) 

DO 180 M=2,MM1 ,2 

AA<1 ,M)»AA<1 ,M>+<ER*EYM<N,M)-EI*EYP<N I M>> 
AA<2,M)=AA<2,M)+<EI*EYM<N,M>+ER*EYP<N,M>> 
AA<3,M)=AA<3,M)+<ER*EZM(N,M>-EI*EZP<N,M>> 
AA< 4 ,M)=AA< 4 ,M> + < El *EZM<N ,M> +ER*EZP<N ,M> ) 
180 CONTINUE 

DO 190 M=1 ,MM,2 

AA< 5, M)=AAC 5 ,M> + < ER*EYM<N ,M> -El *EYP<N ,M) ) 
AA<<S,M)=AA<6,M) + <EI*EYM<N,M) + ER*EYP<N,M>> 
AA<7,M)*AA<7 > M)+<ER#EZM(N,M)-EI*EZP<N f H)) 
AA(8,M)«AA<8|M)+<EI #EZM<N,M>+ER*EZP<N,M) ) 
190 CONTINUE 
200 CONTINUE 
210 DO 280 J=1,NUMPH 
PHI=P1+<J-1)*P3 
PH*=PHI *DTOR 
SP=SIN< PH) 

CP=COS<PH) 

AKY»AKS*SP 

AKYSQ=AKY*AKY 

AKYZ SGNAKYSQ-AKZ SQ 

IF<ABS(AKYZSQ).LT.l .E-05> 80 TO 220 

CC=CO/AKYZSQ 

SAY=1 .0 

IF<ABS<AKY> .GT.l .E-05)SAY=SIN<AKY)/AKY 
CAKY=COS<AKY> 

FA=CC*< CAKZ*SAY-CAKY*SAZ ) 

FB as CC*< SAZ-SAY) 

GO TO 240 

220 IF<ABS(AKY) .LT.0.005) GO TO 230 
IF(ABS<AKZ) .LT. 0.005) GO TO 230 
CC=-2 . 0 *DDSQ/ < AKY **3 > 

SAKY=SIN<AKY) 

CAKY=COS< AKY) 

FA=CC* < SAKY *CAKY -AKY ) 

FB=CC* < AKY*CAKY-SAKY ) 

GO TO 240 
230 FA=C1 
FB=C2 

240 CONTINUE 
FYR=0 .0 
FYI=0 .0 
FZR»0.0 
FZ 1=0.0 



» 
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DO 250 M=1 ,MM 
EM=(M-MT>*AKY 
ER=COS< EM) 

EI=SIN< EM) 

FYR=FYR+<FA*<ER*AA< 1 f M)-EI*AA<2,M> ) 

6c +FB*<ER*AA<5,M>-EI*AA<6,M>)) 

FYI=FYI+<FA#<EI#AA< 1 ,M)+ER#AA<2,M) ) 

6c + FB#< El *AA< 5 ,M)+ER*AA<6 ,M> > ) 

FZR=FZR+<FA*<ER*AA<3,M)-EI*AA<4,M>) 

6c +FB*<ER*AA<7 I M)-EI#AA<8,M>>> 

FZI=FZI+<FA*<EI*AA<3,M)+ER«AA<4,M>> 

6c +FB*< El «AA< 7 ,M) + ER*AA< 8 ,M) ) ) 

250 CONTINUE 

ETR=FZR*CP 
ETI=FZ I *CP 

EPR=FYR*ST+FZR*CT#SP 
EP I =*FY I *ST+ FZ I *CT*SP 
UTH=ETR*ETR+ETI *ETI 
UPB=EPR*EPR+EPI *EPI 
6TH=UTH/U0 
GPH=UPH/UO 
260 DBT=- 150.0 
DBP=-150.0 

I F< GTH . GT . 1 . E-l 5>DBT=1 0 . *AL0G1 0 < GTH) 

I F ( GPH . GT . 1 . E-l 5) DBP=1 0 . *AL0G1 0 < GPH) 
WRITE<12,270)THETA,PHI , DBT , DBP 
270 FORMAT < 4F1 0 .4) 

280 CONTINUE 
290 CONTINUE 
RETURN 

300 WRITE< 12,310) 

310 FORMAT < 1X42HNUMBER OF SAMPLE POINTS EXCEEDS ARRAY SIZE/) 

STOP 
END 
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